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ABSTRACT 

A pre-flight analysis of the directional solidification of BiSn with MEPHISTO-4 is presented, 
simplified Bridgman growth under microgravity conditions is simulated using a two dimensional 
finite element model. This numerical model is a single domain, pseudo-steady state model, and 
includes the effects of both thermal and solutal convection. The results show that for all 
orientations of the applied steady state gravity vector, of magnitude 1 p.g, the directional 
solidification process remains diffusion controlled. The maximum convective velocity was found 
to be 4.424 x 1 O' 5 cm/s for the horizontal Bridgman growth configuration. This value is an order 
of magnitude lower than the growth velocity. The maximum and minimum values of solute 
concentration in the liquid at the crystal-melt interface were 13.867 at.% and 13.722 at.%, 
respectively. This gives a radial segregation value of 1; = 1.046% at the interface. A secondary 
objective of this work was to compare the results obtained to those that consider thermal 
convection only (no solutal convection). It was found that the convective flow patterns in 
simulations which included solutal convection were significantly different from those which 
ignored solutal convection. The level of radial segregation predicted by the current simulations 
is an order of magnitude lower than that found in simulations which ignore solutal convection. 



The final aim was to investigate the effect of g-jitter on the crystal growth process. A simulation 
was performed to calculate the system response to a 1 second, 100 pg gravity impulse acting 
normal to the direction of growth. This pulse is consistent with that induced by Orbiter thruster 
firings. The results obtained indicate that such a gravity pulse causes an increase in the level of 
radial solute segregation at the interface from the steady state values. The maximum value of 
solute concentration in the liquid was found to be 13.888 at.%, the minimum value calculated 
was 13.706 at.%, yielding a radial segregation value of £ = 1.31% at the interface. These values 
occurred 126 seconds after the pulse terminated. Thus it is anticipated that the process will 
remain diffusion controlled even when subjected to such g-jitter. 

1. INTRODUCTION 

The synthesis of advanced materials for the automotive, aerospace, electronic and biomedical 
fields demands high-quality crystals. The compositional uniformity (and hence the quality) of 
such crystals can be profoundly influenced by transport phenomena which occur in the melt 
region. The primary transport mechanism causing these deleterious effects is buoyancy induced 
natural convection in the melt region. The low gravity environment of space offers an 
opportunity to suppress the occurrence of buoyancy induced natural convection. Hence there is 
a great deal of interest in the study of directional solidification of crystals in space. 

The MEPHISTO project [1] is a collaborative United States, French and Australian program of 
space experiments aimed at understanding the fundamental processes involved in crystal growth. 
The space-bome experimental apparatus is a Bridgman type furnace with an isothermal hot zone, 
an isothermal chill zone, and an insulated gradient zone. The apparatus was designed and 
manufactured by the French team. Experiments are conducted on board the Mission-Particular 
Equipment Support Structure (MPESS) located in the Shuttle’s cargo bay. Data acquisition and 
experiment control is performed remotely from the NASA Marshall Space Flight Center. The 
MEPHISTO furnace contains three ingots of BiSn (U.S. experiments) or SnBi (French 
experiments) binary alloy inside fused silica ampoules with a 6 mm inner diameter and a 10 mm 
outer diameter. All three samples are solidified simultaneously. An identical thermal condition 
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is imposed on all three samples by two furnaces. A fixed furnace causes the formation of a fixed 
solid/liquid interface, whilst the moving furnace can facilitate either melting or solidification in 
the sample. This dual furnace arrangement then causes a region of molten alloy which is 
bounded by a stationary solid/liquid interface at one end and a moving liquid/solid interface at 
the other. The apparatus employs the Seebeck technique to determine the moving interface 
temperature in one sample. The fixed interface is required as a reference for the Seebeck 
measurements. The moving interface position and shape is determined in the second sample via 
Peltier pulsing. The final sample is fitted with a resistance measurement circuit which measures 
the interface velocity, and a rapid quenching capability which enables the solute profile in the 
liquid near the interface to be captured for post-flight analysis. After flight, the furnace is 
shipped to France where the samples are extracted. For the U.S. experiments, the samples are 
returned to the U.S. for post-flight analysis by teams at the University of Florida and the NASA 
Lewis Research Center. 

Three MEPHISTO space experiments have taken place previously, on board United States 
Microgravity Payload (USMP) -1, -2, and -3 in October 1992, March 1994 and February 1996, 
respectively. MEPHISTO-4 is scheduled to fly on USMP-4 in November 1997. MEPHISTO- 1 
and -3 were French experiments which investigated the non-faceted solidification of tin doped 
with bismuth solute (SnBi). The MEPHISTO-2 and -4 experiments are U.S. experiments that 
examine the faceted solidification of bismuth doped with a tin solute (BiSn). Alloys based on 
bismuth are selected for their favorable Peltier and Seebeck coefficients [2]. 

The principal aims of the MEPHISTO-2 experiments were (i) to examine the interfacial 
morphological stability threshold of the faceted solid/liquid interface in the absence of 
convection; (ii) to evaluate the effects of interface kinetics on stability during diffusive, fast 
diffusive and convective transport regimes; (iii) to examine Bi growth kinetics as a function of 
supercooling and (iv) to investigate phenomena associated with kinetic roughening [1, 2]. The 
MEPHISTO-2 experiments used a Bi-0.1 at.% Sn alloy. Two successful melting and 
solidification cycles were performed. Subsequent to this, experimental anomalies occurred with 
Peltier pulsing and Seebeck signals. However, fifty-five further melting and solidification cycles 
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were performed. Finally, 15 cm samples were solidified, with different growth rates, for post- 
flight microstructural analysis. The loss of reliable Seebeck and Peltier data compromised many 
of the objectives of the experiment. Despite this, the experimental data resulted in a greater 
understanding of the dominating role of interface kinetics in morphological stability. Stability 
phenomena were observed that had not been previously predicted by theory or measured in 
terrestrial experiments [2], 

The MEPHISTO-4 experiment will build on the results of the MEPHISTO-2 experiment and 
allow information lost due to the anomalies in MEPHISTO-2 to be regained. This experiment 
will use a Bi-1.0 at.% Sn alloy. Modifications to the experimental apparatus include the addition 
to each ampoule of a 2mm inner diameter capillary tube, which extends half-way through the 
process zone. MEPHISTO-2 results indicate that the 6mm cross-section contained a maximum 
of three grains across the diameter. It is anticipated that the 2mm diameter capillary will 
facilitate the growth of a single crystal. Thus the Seebeck data for growth inside the capillary 
can be compared with data obtained for growth beyond the capillary. In this way a comparison 
of the interfacial kinetics of a single crystal can be compared to the polycrystalline case. Similar 
orientation-dependent comparisons can be made using the Peltier sample, while the capillary is 
included in the quenching sample for consistency, and tracking of the resistance change. 

2. MEPHISTO: NUMERICAL WORK 

The MEPHISTO project includes a program of computational modelling of the crystal growth 
process. In particular, the role of convection in the melt region is to be investigated by this 
approach. Accurate experimental determination of convection in metallic melts is very difficult 
to achieve. The opacity and chemical reactivity of the metallic melts make non-intrusive 
measurements difficult. Thus numerical investigation of this phenomena is crucial to the 
complete understanding of the process. Furthermore, the computational models themselves are 
to be improved by a process involving prediction of, and comparison with, the experimental 
results. The aim of this procedure is to develop effective fully three dimensional computer 
simulations of fluid flow related effects [2]. Previously, three-dimensional solutions for 
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Bridgman growth have been limited to steady state growth in succinonitrile (a widely used 
transparent phase change material with properties analogous to metallic materials) [3]. 

One of the critical functions of the supporting computational work is to predict the possible 
effects of convection in the microgravity environment of space. The microgravity levels which 
are imposed on the apparatus occur as a result of tidal forces and Shuttle aerodynamic drag. 
Such accelerations are considered to be quasisteady [4]. The exact orientation of this quasisteady 
g vector cannot be estimated until late in the flight schedule. The effect of transient g impulses 
(caused by crew activity or orbital maneuvering - sometimes called “g-jitter”) need to be 
investigated. 

During MEPHISTO-1, g-jitter as a result of use of the Orbiter Maneuvering Thrusters (OMS) 
caused convection in the melt which was observed in the Seebeck measurements. The effects 
on solute transport of this g-jitter have been modelled by means of a pseudo-steady state model 
[5]. The results exhibited reasonable agreement with the Seebeck data. A more detailed transient 
model for MEPHISTO-1 and -2 has been developed, using FIDAP and including front tracking 
[7]. This model included fully coupled heat transport, solute segregation, fluid motion and phase 
change. A two-step solution approach was needed to implement front-tracking in addition to 
solute segregation. Axisymmetric Pseudo-Steady State Models for MEPHISTO-1 and -2 were 
also developed [7], 

There have been other pre-flight computational analyses in support of MEPHISTO-4. For 
example, sample rehomogenization has been investigated using a variety of finite-difference 
techniques [8]. This was not a simulation of the solidification process, rather the 2D convection- 
diffusion equation was solved for the fluid domain in order to determine the time for solute 
redistribution after a remelting procedure. Both thermal and solutal convection was considered. 
Results indicated that no additional diffusion time was required for the system to reach an 
acceptable level of solute uniformity. By means of a transient, 2D FIDAP finite-element model, 
convection effects at microgravity levels, and in particular the influence of the presence of the 
growth capillary, for MEPHISTO-4 have also been modelled [9]. The model used a fixed-grid 
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approach with the enthalpy method being employed to model the phase change. Temporal 
averaging was used for the apparent heat capacity in the discretized equations. Due to 
computational difficulties introduced by the small partition coefficient for BiSn, the effects of 
solute on density, and thus solutally driven convection, were ignored. 

A scaling analysis of the convection levels by de Groh and Nelson [4] indicate that solutal 
convection effects on solute segregation may be significant. However, thus far, it has not been 
possible to include solutal convection, or even a passive solute into numerical simulations 
involving phase change for MEPHISTO-4. This is due to the difficulties with convergence with 
front-tracking methods and difficulties imposed by the low value of partition coefficient for BiSn 
alloy [7, 9]. The computational modelling presented in this document is intended to examine the 
effects of thermo-solutal natural convection on the MEPHISTO-4 experiments. A simplified 
single domain model will be used. These simulations are intended to supplement ongoing work 
which involves more complex transient, front-tracking models of the process that do not include 
solutal convection. The simplified model employed is a pseudo-steady state model [5, 6], which 
considers only the melt region of the solidification process. Effects of solute rejection due to 
solidification are modeled by an “interface” boundary condition (refer to figure 2) in lieu of 
including the presence of a solid phase. The key assumptions of the model are that the interface 
is planar with its temperature dependent on composition, the diffusion-controlled directional 
solidification is at “steady state," end effects due to the finite length of the ampoule are 
considered negligible, and the furnace translations rate is equal to the crystal growth velocity. 
The constitutive equations for mass, heat, species and momentum transport in the fluid region 
are four coupled, second order, partial differential equations. Solutions to these equations will 
be obtained using the commercial finite element CFD package FIDAP. 

Simulations will be performed for cases involving a steady g vector (representing tidal and drag 
forces) of constant magnitude. The direction at which the g vector acts will be varied for each 
case. In addition, a transient simulation will be performed to investigate the system response 
to g-jitter. A g-pulse will be used which is consistent with that imposed by the thrusters of the 
shuttle during standard orbital maneuvering. 
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3. MATHEMATICAL FORMULATION 


3.1 Governing Equations 

For this work, we shall consider the melt region during the directional crystal growth of a binary 
alloy by the Bridgman process (refer to figure 1). The melt is considered to be a viscous, heat- 
conducting Newtonian fluid subject to thermo-solutal convection. Thermophysical properties 
are dependent on temperature while density variations in the fluid are considered subject to the 
Boussinesq approximation. 

The principle of conservation of momentum is then governed by the Navier-Stokes equation 


( du "1 

— +u-Vu 

> a 


= P 0 8 


i-P t (t-t)-p c (c-c) 


- VP + V 


jp^u+(Vu) n 


( 1 ) 


where u is the velocity vector, p 0 is the fluid density at the reference temperature T 0 and initial 
solute concentration C 0 , T is temperature, g is the gravitational acceleration vector, p T and P c 
the volumetric expansion coefficients due to temperature and solute 1 concentration, p is the first 
coefficient of viscosity and P is the pressure. Since density variations are limited to the 
buoyancy term in equation (1), the governing equation for the principle of conservation of mass 
reduces to the divergence free condition 


Vu = 0 


(2) 


The equation for the transport of energy is a second order partial differential equation 


1 The units of P c , and (C - C 0 ) in P c (C-C 0 ) of (1) were (atomic fraction)' 1 and atomic 
fraction respectively; this was an error resulting in a larger P c (C-C 0 ) product. The proper 
units of P c and C are (volume fraction)' 1 and volume fraction. See Postscript notes, section 7, 
for a discussion. 


7 




in which k is the thermal conductivity and c p is the specific heat capacity. The governing 
equation for the conservation of solute concentration is 

dC 

— +u-VC = V- DVC (4) 

dt ' ' 

where C is the atomic fraction of solute and D is the diffusion coefficient of the solute (Sn) in 
the solvent (Bi). 

The model under consideration is a pseudo-steady state model similar to that by Alexander et. 
al. [6], This is a single domain model which only considers the melt region of the solidification 
process. Thus constitutive equations (1) to (4) are continuous over the entire simulation domain 
and no solid-liquid interface matching conditions need to be developed. However, the 
temperature of the pseudo interface is allowed to change due to changes in solute concentration 
in compliance with the equilibrium phase diagram. 

Now that the constitutive equations have been introduced, initial conditions and boundary 
conditions need to be specified so that the problem becomes fully defined and a solution may be 
sought. 

3.2 Initial and Boundary Conditions 

The pseudo-steady state model (PSSM) is a widely used model for the simulation of steady-state 
Bridgman crystal growth. There are many variants of this model, the simplest being single- 
domain solutions that consider the solid/liquid interface to be flat with a constant melting 
temperature [5, 6]. More complex variants may include the presence of a solid [7]. For the 
present work, we will improve the single-domain steady state model of Alexander [5] and 
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Alexander et. al. [6] by including the effects of a concentration dependent solid/liquid interface 
temperature and solutal convection. This new idealized model is shown in figure 2. The 
assumptions that underpin this new model are 1) that the interface is considered to be planar and 
its temperature concentration dependent, 2) the diffusion controlled directional solidification is 
at “steady state,” making the average composition at the s/1 interface C 0 , 3) end effects due to 
the finite length of the ampoule are negligible, 4) the furnace translation rate is equal to the 
crystal growth velocity, 2 and 5) the (steady state) concentration profile used is a reasonable 
approximation of the profile expected during the (transient) process. Translation of the ampoule 
relative to the furnace is modelled by translating the coordinate axes at a rate equal to the growth 
velocity. This results in a continuous slug of melt moving opposite to growth velocity (-u g ) and 
at the bulk solute composition C 0 being supplied through the inlet. The “inlet” or “far field” 
boundary conditions then are (from ref. [5, 6]) 


at x = x f : u = -u , u = 0, 

x g y 


T = T , 

h’ 



(5) 


in which T h is the hot zone temperature. The effect of solute-poor material solidifying at the 
crystal-melt interface is modelled by the boundary conditions applied at the “interface” boundary. 
Since we have assumed that the crystal growth rate is equal to the furnace translation rate, the 
interface boundary is located at a fixed distance L from the inlet. Since the freezing temperature 
of a binary alloy is dependent on the solute concentration (refer to figure 4), the interface 
temperature will be considered variable with concentration. Thus, the boundary conditions 
applied at the interface boundary become (from ref. [5, 6]) 

at x = X:. : u = - u , u=0, T=T 

X g y 


+ mC , 


dC 

dx 


u 


^(l-k)C 


( 6 ) 


2 Because of the solute build up at the interface, this may not be true during the initial 
transient. However, because of the large temperature gradients, the interface quickly catches 
up with the furnace translation. 
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in which T m is the melting temperature of pure bismuth, m is the slope of the liquidus line (refer 
to figure 4) and k is the partition (or segregation) coefficient for a BiSn alloy. The ampoule 
walls are considered to be impermeable to solute and translate at the growth velocity. 


at y = ±R: u = -u , u =0, 

J X % y 



(7) 


in which R is the ampoule inner radius (half simulation domain height). For the thermal 
boundary conditions, there is an initial adiabatic length near the interface and the remainder of 
the walls are at the furnace temperature 


at y = ±R : 


=n for x s x < x , T = T for x s x s x 

dy i » h . l 


( 8 ) 


in which Xj, x a and x, are the interface, end of adiabatic zone and inlet x locations, respectively. 


The initial conditions employed for the steady state solutions were 
u = -u , u =0, T=T C = C (9) 

x g y h 0 


3.3 Concentration Modification 

Boundary conditions for concentration (5) and (6) implicitly assume that, for concentration, a one 
dime nsional “steady-state” growth mode exists. During steady state growth, the concentration 
of die solid at the interface is equal to the far-field value of concentration (C 0 = 1 at.%). The 
concentration in the liquid is given by the equation (from ref. [10, 11]) 
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( 10 ) 


c =c 

l o 


1 + 


1-k 


expl -u x'/D 

t 


in which x' is the distance from the front into the melt. Note that this analytical solution has 
been obtained by assuming a constant value for diffusion coefficient, D. Due to the low value 
of partition coefficient, k, for the BiSn alloy system considered in this work, the diffusion- 
controlled growth is in the “initial transient” regime throughout the entire experiment. Thus the 
steady state concentration profile in the liquid is never encountered in the experiment. The 
concentration values calculated by the steady state model were found to be unrealistic; in 
particular it predicted concentration values which would lead to the formation of the second 
phase (Sn-rich) at the interface. The model was altered so that it would furnish a more realistic 
concentration profile. The most effective way to achieve this end was to retain boundary 
conditions (5) and (6), and simply alter the far-field concentration value, C 0 , and thus the average 
solid composition at the interface, changing them to C s *. The value of C s * was determined using 
the equation for the concentration at the solid interface, during the initial transient, with diffusion 
controlled growth (from ref. [10, 11]) 


C* 


c 

o 


1 -(l -k)exp[-ku z/D 

v 1 \ « 


(for k < 0.1) (1 1) 


In which z is the distance over which solidification has occurred. The result after growth has 
progressed for 1.5 cm is C s * « 0.4 at.%. The value of C 0 selected was then 0.4 at.%. Figure 5 
is a plot of the steady state profile yielded with this new value. A constant diffusion coefficient 
of D = 3.0 x 10' 5 cm 2 /s (which is the value for D at 550 K, refer to table Al) was used. This 
altered model then constitutes the concentration modified pseudo steady state model (CM-PSSM). 
Its utility is that it enables us to investigate the effects of thermo-solutal convection on the 
solidification process. This compromise was necessary since more realistic models, involving 
a full transient analysis with front-tracking, have not been able to include the effects of solute 
[9], yet it is strongly desired to gauge the effect of solutal gradients on convection, since those 
effects are expected to be of the same order of magnitude as those thermal driven [4,8]. 
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In summary, a new steady-state model known at the concentration modified pseudo steady state 
model (CM-PSSM) has been proposed. The key assumptions employed in formulating this model 
are 1) the interface is planar and its temperature is concentration dependent, 2) a "steady state" 
mode of directional solidification exists, 3) end effects are negligible, 4) the furnace translation 
rate is equal to the growth velocity, and 5) the (steady state) concentration profile used is a 
reasonable approximation of the profile expected during the (transient) process. Note that no 
constitutional supercooling in the melt is incorporated in the model. 

Now that the constitutive equations, initial conditions and boundary conditions have been 
specified for the model, an appropriate numerical solution scheme needs to be determined so that 
a solution may be obtained to the problem posed. 

4. NUMERICAL ANALYSIS 

4.1 Discretization and Solution Scheme 

Numerical solutions are obtained via the commercial finite element CFD code FIDAP (release 
7.6). Four node quadrilateral finite elements were used. Bilinear shape functions for velocity, 
temperature and concentration were employed. Linear interpolation for pressure (discontinuous 
at the element boundaries) was used. Further details of this finite-element formulation may be 
found in [12]. The FEM discretization of equations (1) to (4) results in a set of strongly-coupled 
nonlinear algebraic equations. For the steady g simulations (presented in section 5.1), a Picard 
iteration scheme was used to solve these equations. This involved applying a single iteration of 
successive substitution (large radius of convergence, slowly convergent) followed by nonlinear 
Newton-Raphson iterations (small radius of convergence, quickly convergent). Convergence of 
the iterative scheme is assessed using the following criteria 


||U P -U P ~ 1 || .3 

<; 10 


u‘ 


P- 1 ! 


( 12 ) 
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(13) 


l|R p ll 

«R°II 


<; 10 


-3 


in which u represents the solution vector (of the field variables velocity, temperature and 
concentration), R is the residual vector, p denotes the iteration number and the L 2 norm is defined 
as (from ref. [13]) 



where n is the length of the vector and j is a dummy index. The tolerance value used (10 3 ) is 
an order of magnitude smaller than the FIDAP default value and has been used in previous work 
[3]. Convergence occurred after no more than 5 iterations for all cases studied. Computations 
were carried out on a Silicon Graphics Power Challenge workstation with four 90 MHz 
processors. The CPU burden was typically 50-58 seconds. For the transient simulations 
involving a g impulse (presented in section 5.2), solution at each time step was achieved via a 
Newton-Raphson iteration scheme. The same convergence criteria were employed. Fully implicit 
time integration was used with a variable time step. The minimum and maximum time steps 
used were 0.05 and 8 seconds, respectively. Computations were performed on the same platform 
as for the steady g simulations. A total of 505 time steps were calculated, for a final finish time 
of 1,140 seconds. The CPU requirement was 5,600 seconds. 
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4.2 Grid Independence Study 


The mesh used is shown in figure 3. There are a total of 60 finite elements in the x-direction 
and 18 finite elements in the y-direction. Half of the total number of finite elements are located 
within the adiabatic zone, with the elements being progressively graded as the crystal-melt 
interface is approached. This is to ensure an efficient placement of elements near the crystal- 
melt interface where solutal gradients are large. To determine the adequacy of this discretization, 
simulations were performed using three different mesh sizes, for a test case with a constant 
gravity vector of 1 pg oriented at 0 = 0° (horizontal Bridgman growth). The flow-field results 
indicated that the 60 * 18 discretization was appropriate, since the variables examined did not 
change by more than 1 % at the next highest level of mesh refinement. The results of this study 
are summarized in Table 1. The field variables examined were 

• | U max | , the maximum value of velocity throughout the flow domain. The value of growth 
velocity is subtracted from the velocity vectors used to determine this calculation (refer 
to section 5.1). The location of this quantity is also given. 

• T mid , the value of temperature in the midpoint of the adiabatic zone (x = 4.77305 cm, y 
= 0 cm). This quantity has been used in the literature to test convergence for natural 
convection in a closed cavity [14, 15]. 

• C mid , the value of concentration at the interface boundary and on the centerline. This 
location was selected in the absence of a “natural” test location. It should be noted that 
the concentration profiles were almost identical for all three meshes (indicating that the 
mesh grading is effective in resolving the diffusion-controlled solute boundary layer). 

The CPU burden required is also given in the table. Note the difficulties encountered in 
achieving convergence at the coarsest mesh. 
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Table 1. Summary of grid independence results (quantities in brackets indicate the percent 
change between the value found using the present mesh and at the next most refined mesh). 


Mesh 

o 

X 

X 

■B 

£ 

ID 

|U m J location 

T m id 

c 

^mid 

CPU 


(cm/sec) 

(x cm, y cm) 

(K) 

(at. %) 

(sec/# iterations) 

30 x 9 

4.514 

(5.123, 0.1667) 

755.0 

13.7921 

101/29 


(2.03%) 


(0.453%) 

(0.001%) 


60 x 18 

4.424 

(5.132, 0.1667) 

751.6 

13.7919 

58/5 


(0.6827%) 


(0.173%) 

(0.001%) 


120 x 36 

4.394 

(5.137, 0.1667) 

750.3 

13.7918 

468/5 


5. RESULTS 

5.1. Steady Gravity Vector 

Simulations were performed for a g vector of constant magnitude 1 pg at a variety of orientations 
(the definition of angle 0 is provided in figure 2). The orientations considered were 0 = 0° 
(horizontal Bridgman growth), +22.5°, +45°, +67.5°, +90° (vertical Bridgman growth, crystal at 
top) -45° and -90° (vertical Bridgman growth, crystal at bottom) degrees. The most pertinent 
results were for the 0 = 0°, +45°, +90° and -90° configurations because residual gravity in 
MEPHISTO-4 is expected to be 1 pg at +45° and previous MEPHISTO modeling has considered 
gravity vectors acting at 0° and -90°. These will be presented in depth in this section while plots 
of the results for the remaining configurations are provided in appendix A. In addition, 
supplementary simulations were performed which excluded the effect of solutal convection and 
considered thermal convection only. The orientations considered were 0 = 0°, +45°, +90° and 
-90° degrees. The purpose of these supplementary calculations is to provide a comparison with 
results obtained in the cases which included solutal convection, and hence aid in the 
determination of the importance of solutal convection on the interface solute concentration levels. 
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Thermophysical properties and other constants used in obtaining the solution are provided in table 
A1 in the appendix. Note that the x locations of the interface, end of adiabatic zone and inlet 
were taken to be x ( = 3.94805 cm, x a = 5.59805 and x, = 11.25 cm. These values were consistent 
with the simulations for MEPHISTO-4 in ref. [9]. 

Figure 6 (a) is a plot of velocity vectors for the case of 0 = 0’ (horizontal Bridgman growth). 
Note that the growth velocity component has been subtracted from these plots for clarity. Also 
note that the entire domain is not shown; the x range of the plots is smaller than the cavity length 
in order to highlight details near the interface. The main feature of the results is that there are 
two counter-rotating convective cells in the flow domain. The primary convective cell, which 
is driven by thermal gradients, rotates in a counter-clockwise fashion. The length of this 
convective cell is approximately 2.1 cm resulting in aspect ratio (height: length) of 0.29. The 
maximum velocity in the domain occurs in this thermal convective cell. The magnitude of the 
maximum velocity is 4.424 x 10 5 cm/s (an order of magnitude lower than the growth velocity) 
at a location of (5.132, 0.1667) cm. As will be shown later, at these velocities diffusive transport 
is dominant over convective transport. The secondary convective cell, driven by solutal 
gradients, is near the interface. This cell is 0.172 cm long with an aspect ratio (height: length) 
of 3.49. The velocities in this cell are smaller than in the primary cell. In particular, the 
velocities are very small near the interface, for example, the velocity at (x = 3.9462 cm, y = 0 
cm) is only 3.873 x 10' 6 cm/s. Figure 6 (b) is a plot of the density variation in the buoyancy 
term of equation (1). This plot makes the net driving force behind the convective motion seen 
in figure 6 (a) readily apparent. At the interface, solute rejection leads to high levels of solute. 
This opposes the tendency of the fluid density to increase in response to the low temperatures 
near the interface. The net result is that the fluid is below the reference density at the interface 
(note that Sn is lighter than Bi). Moving away from the interface, the solute concentration 
decays while temperatures steadily increase. The net result is a region of maximum density 
located between the two convective cells seen in figure 6 (a). Further away from the interface, 
the density decreases linearly in the adiabatic zone. Near the end of the adiabatic zone the 
density contours exhibit curvature. Past the adiabatic zone the density takes on a minimum 
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value. Figure 7 (a) is a plot of solute concentration contours for this gravity orientation. Note 
the finer scale used on the x-axis (x range to 6 cm). The main feature to note is that the solute 
concentration lines remain flat. This indicates that convective effects on the solute boundary 
layer are negligible and so the growth is dominated by diffusion. Figure 7 (b) is a plot of the 
temperature contours for the same case. Again, convection has not significantly distorted the 
isotherms. Note also the curvature of the isotherms at the end of the adiabatic zone (x = 5.59805 
cm). These are due to heat being conducted in from the hot zone. The minimum temperature 
at the interface is 512 K; this is lower than the reference value of 544.3 K since the solute 
concentration at the interface has suppressed the melting temperature. The interplay of 
concentration levels indicated in figure 7 (a) and the temperatures shown in figure 7 (b) result 
in the density contours shown in figure 6 (b). 

Figure 8 (a) is a plot of velocity vectors for a gravity orientation of 0 = +45”, i.e. equal 
components of gravity acting in the negative y direction and horizontally away from the interface 
in the positive x direction. As for the horizontal growth case detailed above, the most important 
feature is the presence of two counter-rotating convective cells. The character and dimensions 
of the cells are identical to those in the horizontal case, however the velocities are smaller in 
magnitude. The maximum velocity occurs in the primary convective cell. Its magnitude is 3.132 
x 10' 5 cm/s at a location of (x = 5.132 cm, y = 0.1667 cm). The ratio of the maximum velocities 
for this orientation to that found for horizontal Bridgman growth is 0.707. This is also the ratio 
of the y gravity components, highlighting the direct relationship between y gravity component 
and the resulting level of convection. There is no discernible difference between the plots of 
density, solute concentration and temperature contours for this case and those provided in figures 
6 (b), 7 (a) and 7 (b) and so they have been omitted. Figure 8 (b) is a plot of velocity vectors 
for a gravity orientation of 0 = +90”, i.e. anti-vertical Bridgman growth with the gravity vector 
pointing away from the crystal and into the melt. For this orientation, concentration gradients 
are stable while temperature gradients are unstable. As there is no longer a y gravity component, 
the character of the convective motion has changed compared to the convective motions exhibited 
in the previous two cases. For the 0 = +90“ configuration, two convective cells exist in the 
region at the end of the adiabatic zone. The cells are mirror images of each other about the x- 
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axis. The velocities at the x-axis (y = 0) are in the positive x direction. The maximum velocity 
is 3.0498 x 10 -6 cm/s at a location of (x = 5.718 cm, y = 0 cm). This is two orders of magnitude 
lower than the value of growth velocity, and approximately 12 times slower than the maximum 
velocity for the 0 = 0° case. Because of the low gravity magnitude no Rayleigh-Benard style 
convective instabilities [16] are observed. Since there are no thermal or solutal gradients normal 
to the direction of the gravity vector, convection is absent near the interface. Thus conditions 
at the interface are completely diffusion controlled. Again, there is no discernible difference 
between the plots of density, solute concentration and temperature contours for this case and 
those for 0 = 0° so they have been omitted. Figure 8 (c) is a plot of velocity vectors for a 
gravity orientation of 0 = -90”, i.e. vertical Bridgman growth with the gravity vector pointing 
down at the crystal. For this orientation, concentration gradients are unstable while temperature 
gradients are stable. The important detail to note is that the convective motions operate in the 
reverse sense to those for the 0 = +90' case, such that the flow at the centerline is in the negative 
x direction. 

The next step in our analysis is to assess the impact of these solutal-convective motions on the 
solute segregation at the interface (and hence the solute segregation in the solid). Figure 9 shows 
solute concentration values across the interface (x = x, = 3.94805 cm), for each gravity angle 0 
tested. For the anti-vertical and vertical Bridgman growth configurations (0 = ±90°), 
concentration remains constant since conditions at the interface are completely diffusion 
controlled. For other values of 0, the presence of convective transport causes concentration to 
vary slightly across the interface. The minimum and maximum values of concentration occur 
at y = -0.3 cm and y = +0.3 cm, respectively. As 0 moves from +90° to 0°, the difference 
between maximum and minimum concentrations increases. The maximum disparity between 
minimum and maximum concentrations occurs for the 0 = 0' case, corresponding to the 
maximum level of convection. For this orientation, C max = 13.867 at.%, C min = 13.722 at.%. 
These values correspond to 0.402 at.% and 0.398 at.% in the solid. 
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Figure 10 is a plot of radial solute segregation, £, against gravity angle, 0. The radial solute 
segregation is defined as 
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C av is taken to be the average value across the interface (a lateral average [7]). The maximum 
value of segregation is 1.046% for 0 = 0”. The authors consider this value of segregation to be 
negligible. This confirms that the levels of convection are small, the process is diffusion- 
dominated. As angle 0 increases, the segregation level decreases. Note that the level of 
segregation decreases much more rapidly with 0 as 0 - +90”. At 0 = +90”, the pure diffusion 
orientation, there is zero segregation. 

What remains now is to compare and contrast these results to a case in which solutal convection 
is ignored. Figure 1 1 (a) is a plot of velocity vectors for a simulation with 0 = 0” and thermal 
convection only. In contrast to the results for thermo-solutal convection, there is a single, 
dominant counter-clockwise rotating convective cell. The length of this convective cell is about 
2.25 cm with an aspect ratio (height: length) of 0.266. The maximum velocity is 4.427 x 10' 5 
cm/s (slightly larger than the thermo-solutal case) acting at a location of (x = 5.1324 cm, y = 
0.1667 cm), which is identical to that found in the thermo-solutal case. The maximum velocities 
within the solutal boundary layer are 5.450 x 10' 6 cm/s and 1.969 x 10 5 cm/s for the cases 
including and excluding solutal buoyancy, respectively. Figure 11 (b) is a contour plot of the 
density in the buoyancy term of equation (1). Here, there is no contribution from solute, so the 
density varies linearly with temperature throughout the adiabatic zone to the interface. Figure 
12 is a trace of density values along the centerline extracted from figures 6 (b) and 1 1 (b). Note 
the profound influence of solute on the density near the interface, when compared to the thermal 
convection only case. It is these density gradients which are the cause of the convective patterns 
observed in each case. Figure 13 (a) is a plot of concentration contours for the thermal 
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convection only case. The important detail to note is that even though |U max | is about the same 
for the cases including and excluding solutal convection, the maximum velocity within the solutal 
boundary layer is 3.61 times greater in the thermal convection only case; the higher level of 
convection has resulted in the concentration contours being slightly distorted. The contours now 
take on a very gentle "reverse-s” shape. Figure 13 (b) is a plot of isotherms for this case. Note 
that the minimum temperature (510 K) is lower than that found from the thermo-solutal case (512 
K). ' This is expected and is due to the greater amount of solute transport by convection at the 
interface (recall that the melting temperature at the interface is dependent on the value of 
concentration). The isotherms near the interface are still flat. Figure 14 is a plot of 
concentration traces across the interface for the thermo-solutal convection case and the thermal 
convection only case. A diffusion-only case (0 = ±90") is included to provide a basis for 
comparisons. Note the following dramatic difference between the concentration traces predicted 
by the two simulations. First, the locations of the maximum and minimum values of 
concentration have been transposed, such that the segregation is in the opposite sense. Second, 
the higher levels of convection near the interface have resulted in a higher maximum 
concentration (14.51 vs. 13.87 at.%) and a lower minimum (13.02 vs. 13.72 at.%). The value 
of radial solute segregation, £, is an order of magnitude higher, at 11.5%; and is considered large 
enough to represent a departure from diffusion dominated growth. 

5.2 “g-jitter” Analysis 

Transient simulations (including both thermal and solutal buoyancy) were performed to examine 
the system response to the application of a single g impulse. The g impulse vector was at an 
orientation of 0 = 0°, with a magnitude of 100 pg and a duration of 1 second. This corresponds 
with standard orbital maneuvering. In addition, a steady g vector, at an orientation 0 = 0° and 
a magnitude of 1 pg, was applied for the entire duration of the simulation. This corresponds to 
the quasisteady g effect of tidal and aerodynamic drag forces as considered in section 5.1. The 
initial condition used for the solution was the steady state solution obtained for the corresponding 
case in section 5.1. The g impulse was applied from a simulation time of t = 0.0 to t = 1.0 
second. 
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Figure 15 is a plot of velocity vectors at four selected times throughout the simulation. Figure 
15 (a) is at a simulation time of 0.5 seconds, i.e. half way through the gravity pulse. The level 
of convection has increased significantly over the steady state result. The value of maximum 
velocity is 4.080 * 10' 4 cm/s, which is higher than the growth velocity. The location of this 
maximum velocity has shifted to (x = 4.861 cm, y = 0.2333 cm). This is further from the 
centerline and closer to the interface than for the steady state case. Indeed, the entire convective 
motion is much stronger, as can be gleaned by perusing the plot. This trend continues for the 
plot shown in figure 15 (b). This figure is the velocity field at a time of t = 1.0 second, i.e. just 
as the pulse terminates. This is the time at which maximum velocities are found. The value of 
maximum velocity is 6.309 * 10"* cm/s, which is 1.87 times the growth velocity. The location 
at which the maximum velocity acts is (x = 4.945 cm, y = 0.2333 cm). Figure 15 (c) shows the 
velocity vectors after 25.2 seconds. At this time, the strength of convection has decreased 
appreciably - the maximum value is 6.996 x 10' 5 cm/s at (x = 5.237 cm, y = 0.1667 cm). The 
convective motions are similar to the steady state case although the magnitudes are higher. The 
final plot, figure 15 (d) is for a time of 126 seconds. At this time, the convective motion has 
almost been restored to steady state. The maximum velocity is 4.426 x 10‘ 5 cm/s at (x = 5.132 
cm, y = 0.1667 cm). This value is only 0.1% higher than that which occurs at steady state and 
is at the same location. 

Figure 16 is a plot of radial solute segregation at the interface (as defined by equation 15) against 
time. The maximum value of segregation does not correspond to the time at which convection 
is a maximum; rather the maximum rate of increase in segregation occurs when convection levels 
are at their highest. The increased convection levels in the flow field from 0 to 126 seconds lead 
to increasing levels of radial solute segregation. The maximum segregation is found to be 1.31% 
, at a simulation time of t =126 seconds, when the level of convection is restored (almost) to the 
steady state value. After this time the level of segregation slowly decreases as the solute field 
slowly drifts toward its steady state configuration. Even though £ has increased by about 30% 
due to the impulse, this level of segregation is still insignificant and would be within the scatter 
and uncertainty of a physical measurement. Post-flight measurements on a Bi- 0.1 at%Sn ingot 
solidified in space with similar g and g-jitter conditions have shown composition to increase from 
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about 0.01 at.% to 0.1 at.% along the centerline as a result of the g-jitter impulse [17]. Such a 
dramatic composition increase has not been predicted by the present work. Since the simulations 
performed were subject to steady state assumptions (refer to section 3), the authors suspect that 
all the important transport phenomena have not been modelled adequately. To do this a more 
detailed model is required which includes the presence of a solid phase and a moving solid/liquid 
interface. Such a model is unencumbered by steady state assumptions and should furnish results 
with better agreement with the experimental data. 

Figure 17 is a plot of solute concentration traces across the interface at the simulation times 
considered above. The t = 0.5 and t = 1.0 second traces are almost indistinguishable from the 
steady state traces. Even though the convection is much higher at these times, the convective 
fields have not established themselves for a long enough period of time to have an impact of the 
solutal field. The t = 25.2 and t = 126 second traces show deviation from the steady state. Time 
has allowed the increased convection to alter the solute field. Also worthy of note is the 
asymmetry of the traces at the later times. The departure from the steady state case is less in the 
y = -0.3 to 0 cm section of the interface when compared to the y = 0 to 0.3 cm section of the 
interface. This could possibly indicate that the convective cell driven by thermal gradients is 
slightly dominant over the cell driven by solutal gradients at increased levels of convective 
strength. 

6. CONCLUSIONS 

The aim of these simulations was to quantify the effect of thermo-solutal convection on the 
MEPHISTO-4 space experiments. The results show that for all orientations of the steady state 
gravity vector of magnitude 1 pg the directional solidification process remains diffusion 
controlled. The maximum convective velocity was found to be 4.424 * 10' 5 cm/s for the 
horizontal Bridgman growth configuration (0 = 0°). This value is an order of magnitude lower 
than the growth velocity of 3.38 * 10 -4 cm/s. The maximum value of radial solute segregation 
(as defined by equation 15) was found to be i = 1.046%. The radial segregation was such that 
the concentration at the interface decreased in the direction of applied gravity. A secondary 
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objective of this work was to compare results obtained including solutal convection with those 
that consider thermal convection only. This work indicated that the addition of solutal convection 
had a profound effect of the convective flow patterns and the resulting levels of radial 
segregation at the interface. Simulations that ignored solutal convection were found to 
overpredict the level of radial segregation by an order of magnitude and to reverse the 
concentration profile at the interface, with concentration at the interface increasing in the 
direction of applied gravity. 

The final component of the present study was to examine the effect of g-jitter. For a 100 pg 
gravity pulse of duration 1 second and acting normal to the direction of growth, the process 
remains diffusion controlled. The level of radial solute segregation increased to a peak value of 
\ = 1.31%. The maximum level of convection was found to be 6.309 x 10 -4 cm/s at a 
simulation time of t = 1 sec (i.e. at the time when the pulse terminated). After 126 seconds, the 
level of convection had returned to a value of 4.426 * 1 O' 5 cm/s which is within 0.1% of the 
initial steady state value. This was the time at which maximum radial solute segregation was 
observed. After this time, the level of radial segregation slowly decreased, to a value of 1.11% 
at 1140 seconds (the simulation finish time). 

7. POSTSCRIPT 

There is a subtle and significant nuance to the use of P c and C in equation (1). As defined, p c 
is equal to dp/dC divided by p 0 . In practice however, P c is determined as (p^p^/p^Q-CJ 
where pj and p 2 are the densities of the pure solvent and solute respectively, and C t and C 2 are 
the compositions of the pure solvent and solute respectively. Note that C, and C 2 are always 
equal to 0 and 100% respectively. In using (pj-p^/p^Cj-C,), a linear relation is assumed 
between density and composition in the units of C. This assumption is valid, in general, when 
the units of C are in volume percent (or fraction) [18]. We used atomic fraction which 
exaggerated the influence of solute in the density of the liquid in (1). However, since the flow 
is nearly negligible, and growth was determined to be dominated by diffusion, this error is not 
considered significant to the major conclusions of this work. 
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Figure 1. Schematic of Bridgman crystal growth process. 
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Figure 2. Computational domain for pseudo-steady state model, including pertinent boundary 
conditions. Angle 0 is defined as the angle from the y-axis to the gravity vector, counter- 
clockwise positive. This results in the gravity force components g, = g sin 0, g y = -g cos 0. 
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Figure 3. The finite element mesh employed for the analysis. Note the grading of the mesh as it approaches the interface. This is 
to adequately resolve solutal concentration gradients near the interface. The aspect ratio (length:height) of the cavity is 12.17:1. 
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Figure 4. Bi rich portion of the linearized BiSn binary alloy phase diagram. The change in 
melting temperature of the alloy with increasing tin composition is given by the slope of the 
liquidus line, m. The composition of the liquid at the interface (C L *) and the composition of the 
solid at the interface (Q*) is related by the segregation coefficient, k. 
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Figure 6. (a) Velocity vectors for 0 = 0" simulation, including the effects of buoyancy due to 
temperature and solute. The growth velocity, u g has been eliminated from the x-component of 
velocity for this and subsequent velocity plots so the true nature of the convective flow can be 
revealed. Two convective cells can be discerned, one small, clockwise rotating cell next to the 
interface and a large counter-clockwise rotating cell, (b) A contour plot of the value of density 
used in the buoyancy term of equation (1). This indicates the gradients which motivate the 
convective flow seen in (a). 
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Figure 7. (a) Contour plot of solute concentration. Note the high gradients near the interface 
and that the lines have barely deflected from the vertical, (b) Contour plot of temperature 
throughout the simulation domain. Note the suppression of interface temperature from the 
zero concentration value of T m = 544.3 K. 
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Figure 8. Vector plots for (a) 0 = +45 , (b) 0 = +90 and (c) 0 = -90 cases. Note the different 
scales for the velocity vectors. 
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Figure 9. Solute concentration plot across the interface, in the liquid, for each orientation. The 
“bottom” of the interface is to the left (y = -0.3 cm), the “top” is to the right (y = +0.3 cm). 
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Figure 10. Solute segregation vs gravity vector angle. Segregation is at a maximum for 
horizontal Bridgman growth (0 = 0°). 

0 - 0, 1 jig, no solutal conv. 


— -*• 1 cm - 7 2030 x 10* cnV« 

4 5 6 7 8 



DENSITY: 9 53 9.61 9.70 9.78 9.86 9.95 10.03 1 0.11 



4 5 6 7 8 

x(cm) 


Figure 11. (a) Velocity vectors for 0 = 0' simulation, ignoring solutal convection. Only one 
convective cell exists, (b) A contour plot of the value of density used in the buoyancy term 
of equation (1). This indicates the gradients which motivate the convective flow seen in (a). 


33 





4 5 6 7 8 



x(cm) 


Figure 12. Trace of density along the centerline taken from figures 6 (b) and 1 1 (b). The 
presence of volumetric expansion due to solute concentration can clearly be seen. 
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Figure 13, (a) Contour plot of soiutal concentration. The lines exhibit some deflection from the 
vertical, (b) Contour plot of temperature throughout the simulation domain. Note the minimum 
interface temperature is lower (510 K) than that for the case with soiutal convection (512 K). 
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Figure 14. Comparison of the calculated solutal distributions at the interface. The results are 
significantly different for the two cases. 
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Figure 15. Velocity vectors for a simulation time of (a) t = 0.5 sec, (b) t = 1.0 sec, (c) t = 25.2 
sec and (d) t = 126 sec. Note the different scales used for the vector magnitudes. 
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Figure 16. Plot of radial segregation, across the interface agianst simulation time, after a g- 
jitter impulse of 100 [ig for 1 second starting at t = 0. 
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Figure 17 • Solute concentration traces at the interface at various times. Note that the solid line 
for the steady state solution effectively blanks out the lines for t = 0.5 and t - 1.0 second. 
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Figure Al. Velocity vectors for simulations with solutal convection for (a) 0 = +22.5’ and (b) 
0 = +67.5" cases. The character of convective motion is similar in both cases, but note the 
difference in magnitude of velocities. Contour plots of density, temperature and solute 
concentration are not discemably different to those shown in figure 6b and 7, and so are not 
provided. 
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Figure A2. Plot of the maximum magnitude of velocity vs gravity vector angle. Maximum 
convective motion occurs during horizontal (0 = 0°) Bridgman growth. 
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Figure A3. Velocity vectors for simulations without solutal convection for (a) 0 = +45*, (b) 
0 = +90* and (c) 0 = -90* cases. Figures A3 (b) and (c) are very similar to figures 8a and b, and 
the phenomena at the interface are absent from the effects of convection. 
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Figure A4. Concentration contours for simulations without solutal convection for (a) 0 = 
+90' and (b) 0 = +45° orientations. Contours for the -90' case are identical to the +90' case 
and so are not presented. 



Figure A 5. Temperature contours for simulations without solutal convection for (a) 0 = +90° 
and (b) 0 = +45° orientations. Again, contours for the -90' case are superfluous. 
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Figure A6. Solute concentration traces at the interface for simulations without solutal 
convection. Note the sinusoidal shape of the curves. 
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Figure A7. Solute segregation values for each simulation without solutal convection. 
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Table Al. Thermophysical properties and other pertinent physical values. 
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